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DISTRIBUTED PARAMETER ANALYSIS OF PRESSURE AND FLOW 
DISTURBANCES IN ROCKET PROPELLANT FEED SYSTEMS 
by Robert G. Dorsch, Don J. Wood,* and Charlene Lightner 
Lewis Research Center 

SUMMARY 

A digital distributed parameter model for computing the dynamic response of propel- 
lant feed systems is formulated. The analytical approach used is an application of the 
wave- plan method of analyzing unsteady flow. Nonlinear effects are included. The model 
takes into account locally high compliances at the pump inlet and at the injector dome re- 
gion. Examples of the calculated transient and steady -state periodic responses of a sim- 
ple hypothetical propellant feed system to several types of disturbances are presented. 
Included are flow disturbances originating from longitudinal structural motion, gimbaling, 
throttling, and combustion- chamber coupling. The analytical method can be employed 
for analyzing developmental hardware and offers a flexible tool for the calculation of un- 
steady flow in these systems. 


INTRODUCTION 

A variety of analytical models have been employed for calculating the low-frequency 
(0 to 500 cps) dynamic response of rocket propellant feed systems. Early theoretical in- 
vestigations (e. g. , refs. 1 to 3) were primarily concerned with instabilities resulting 
from the coupling of propellant feed system flow perturbations with combustion chamber 
pressure oscillations. More recently instabilities arising from propellant flow perturba- 
tions driven by longitudinal structural oscillations have received attention (refs. 4 to 6) . 

In this type of instability the propellant flow perturbations result in engine thrust pertur- 
bations, which in turn feed back through the structure to close the loop. The authors of 
references 1 to 6 have employed small perturbation techniques and have linearized system 
impedances at mean operating conditions. These simplifications are usually valid for the 
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purposes of making closed-loop stability predictions where the growth or decay of small 
(or incipient) sinusoidal pressure and flow perturbations is considered. 

Because of the generally nonlinear nature of liquid flow systems, linearized models 
cannot be used to calculate the amplitudes and wave shapes of large amplitude pressure 
and flow disturbances of the type that actually occur in propellant feed systems when in- 
stabilities are present. Further, a linearized model cannot be used to calculate the sys- 
tem response to throttling, startup, or shutdown. Nonlinear models have been formu- 
lated to calculate the dynamic and transient responses of rocket engine and propellant 
systems. These models usually employ finite- difference techniques with considerable 
lumping of parameters and usually have limited versatility because they were formulated 
to solve specific problems. 

An analytical investigation was therefore undertaken at the Lewis Research Center to 
develop a general nonlinear distributed parameter model for calculating unsteady flow in 
liquid rocket propellant systems. The objective was to develop a model that can be used 
to calculate both transient and dynamic responses of a liquid rocket propellant feed sys- 
tem to a variety of different types of input disturbances including structural motion, A 
pulse synthesis method fundamentally similar to a method of characteristics solution 
(ref. 7) was chosen for the analytical approach. This method, called the wave plan, was 
developed in the first phase of the investigation and is described in detail in reference 8. 

In the wave-plan method of analyzing liquid filled systems the response of system 
elements or components to incremental pressure pulses is determined. The incremental 
pressure pulses are associated with incremental flow changes, which can originate in a 
fluid system in a variety of ways including mechanical motion of the components. The 
analyses for the response of individual system elements are incorporated into digital com- 
puter subroutines. These subroutines are then combined to form digital computer models 
of specific systems. 

This report illustrates how the wave-plan method can be applied to a liquid rocket 
feed system. Equations for the major feed system components (tank, pump, injector) are 
derived and programmed as digital computer subroutines. These equations and subrou- 
tines along with those of reference 8 are used to analyze a simple monopropellant (or one 
side of a bipropellant) feed system. Examples of the calculated responses of the hypo- 
thetical feed system to several types of disturbances are presented. Included are flow 
disturbances originating from longitudinal structural oscillations, throttling, gimbaling, 
and combustion- chamber coupling. In addition, the effect of local cavitation-induced 
pump inlet compliance on the suction line response is examined. Pump inlet pressure 
and flow and combustion- chamber pressure perturbation wave shapes are calculated and 
presented to illustrate the nonlinear nature of the fluid system. Finally, the effect of 
pump-inlet flow perturbation suppression devices on the feed system response is ex- 
amined. 
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ANALYSIS 


Feed System Model 

A schematic of the idealized system chosen for analysis is shown in figure 1. The 
primary elements or components of the propellant feed system are as follows: 

(1) A propellant tank with provision for longitudinal structural motion input at the 

tank outlet junction with the feed line. The motion input provides for feed system 
coupling when there is mechanical motion of the structure due to longitudinal os- 
cillation of the rocket. The upper boundary of the liquid in the tank is assumed 
to be at constant pressure. 

(2) A constant-diameter propellant feed line (suction line) connecting the supply tank 

and the pump. The motion of the line itself will be ignored in the examples of 
this report although it is recognized that very small disturbances can be gener- 
ated because of viscous effects. (For tapered lines or lines with area changes 
the motion of the line itself must be included. ) 
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(3) A pump with provision for longitudinal structural motion input at the junction with 

the suction line. A device is located at the inlet to the pump in order to simulate 
the cavitation -induced compliance in the inducer region. This device is repre- 
sented by an analytical equation. 

(4) A fix device (perturbation suppressor) located at the pump inlet station. In the 

analysis this device can be mathematically combined with or added to the pump 
cavitation compliance, or it can be treated separately when a mechanical-motion- 
compensating type of fix device (e. g. , constant- volume bellows) is considered 
(ref. 6). 

(5) A pump discharge line. Since the vertical distance from pump to injector is usu- 

ally small, the discharge line will be assumed to be horizontal. 

(6) A nonlinear throttle valve located in the discharge line. 

(7) A nonlinear orifice at the termination of the discharge line. This orifice simu- 

lates the engine injector. The orifice does not move in the horizontal direction. 
The compliance of the injector dome area is simulated by a device (represented 
by an equation) just upstream of the injector. The orifice exit pressure (thrust 
chamber pressure) is a time -dependent function of the mass discharge of fuel 
into the thrust chamber. This time dependence includes both a dead time and an 
exponential delay associated with the gas residence time. 

(8) The viscous resistance of the pump suction and discharge lines simulated (when 

desired) by internal nonlinear (square law) orifices, which approximate turbulent 
flow friction effects. 

Equations describing the response of each of the system elements or components to 
elastic pressure waves resulting from incremental flow changes must be written. These 
equations must be written at the discontinuities in the system (e. g. , the junction of the 
suction line with the pump) because the wave -plan analysis supposes a system composed 
of a discrete number of discontinuities connected by lossless line segments, which serve 
only to transmit pressure pulses. If any of the elements or components are moving with 
respect to the primary reference system lor the elastic waves (a reference system that 
moves with the vehicle but does not oscillate), the motion must be taken into account when 
writing the continuity equations at the discontinuities. In this example specific provision 
is made for including the effects of mechanical motion at both the tank outlet and at the 
pump inlet junctions with the feed line. 


Component and Junction Equations 

Equations for the response of the feed system components to incremental pressure 
pulses will be derived. The general case in which pressure pulses arrive simultaneously 
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from both sides of a junction or dis- 
continuity will be considered (ref. 8). 

Ta nk outlet . - Figure 2 shows 
pressure head and velocity conditions 
before and after impingement of 
elastic- wave pressure pulses at a 
moving tank outlet. 

The tank exit is streamlined to 
minimize viscous losses, for flow 
going into the feed line, and these losses may be neglected. The energy equation for this 
case after wave action is 



Ffgure 2. - CondTIIoris at tank outlet before and after wave action. 


Ml 

2g 


H n = 


_22 

2g 


+ H, 


22 


( 1 ) 


(Symbols are defined in appendix A. ) 

Momentum relations across the wave fronts give the following equations for the pres- 
sure head changes associated with the waves: 


AH, = (V' - V,) 


AH u ’ V <v ' ' v n ) 

Cn 

AH 2 = _£ (V 2 - V’*) 


ah 22 = ~2(v 22 - V") 
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(See fig. 2 for definition of V* and V'*.) Eliminating V' and V" from these equations 
gives 
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Figure 3. - Continuity conditions at tank outlet. 


H 22 “ H 2 + AH 2 + AH 22 ( 5 ) 

Continuity relations are derived by referring to figure 3, which shows the relative 
positions of the tank outlet a short time At apart. The average velocity of the tank out- 
let during this time is V^. 

The net volume flow across section AA minus the net flow out of section BB must 
equal the total storage in time At. (Sections AA and BB move with the vehicle but do 
not oscillate. ) Mathematically, the flow difference is given by 

(V n A t ) At - (V 22 A f ) At = V t (A t - A f ) At 

Eliminating the incremental time and expressing the result in terms of the ratio of 
the feed line to tank area Rj give: 

V n = V t (l - Hi) + y 2 2 n i < 6 > 

The difference in pressure across the tank outlet after the wave action is given by 

H 11 - h 22 ■ H 1 + iH l + iH ll - h 2 ’ iH 2 ‘ AH 22 < 7 > 

Solving equations (1), (2), (3), (6), and (7) simultaneously gives the velocity in the 
discharge line after the wave action in terms of the initial pressure and velocity condi- 
tions, the velocity of the tank outlet, and the magnitude of the impinging waves: 

a V 22 + b V 22 + c = 0 (8) 

where: 
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a = 



b = 


C 1 R 1 

g 



C 2 V 2 

c = H 2 + 2AH 2 - H x - 2AH 1 - — ^ 


C 1 V 1 


g 


' + C l V t 




The positive root of equation (8) gives the propellant velocity in the tank after wave 
action. The propellant velocity in the feed line is given by equation (6). The magnitude 
of the generated pressure waves are given by equations (2) and (3), and the tank and line 
pressures by equations (4) and (5). 

Pump . - Because of the present lack of dynamic data for pumps, the relations between 
net head, suction head, and flow rate for the pump will be assumed to obey the steady - 
state relations at any given instant. Figure 4 shows plots of typical pump characteristics. 

In many cases the relation between the 
net increase in static pressure head across 
the pump and the flow rate Q can be ac- 
curately described by a quadratic expres- 
sion: 




Figure 4. - Typical pump characteristics. 


AHp = A + BQ + CQ* 


(9) 


To obtain the coefficients, A, B, and C the conventional total head rise against flow curve 
(fig. 4) is first converted to static pressure head rise against flow corresponding to the 
particular diameter discharge and suction lines being used. The coefficients. A, B, and C 
may also be expressed as functions of suction pressure H g and rotational speed, which 
makes equation (9) a general expression for static head rise as a function of suction head, 
flow rate, and rotational speed. For values of suction head greater than (fig. 4) and 
fixed rotational speed the coefficients A, B, and C remain essentially constant. 

For the purpose of determining the effects of pressure wave pulses on a moving pump, 
the use of coefficients based on the suction head prior to impingement of the waves suf- 
fices because only incremental pressure pulses are dealt with in the wave -plan method. 

Under most operating conditions a certain amount of cavitation occurs in the inducer 
region of the pump inlet because of high-speed impeller motion in an area of relatively 
low pressure. The presence of this vapor in the liquid results in a locally high fluid com- 
pliance at the pump inlet. This phenomena is very difficult to express in exact mathe- 
matical terms; however, the effect of this lumped compliance must be included in the dy- 
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namic analyses. This effect will be analytically included by attaching a hypothetical 
accumulator at the pump inlet. 

The compliance a' of the hypothetical accumulator is defined by 

dH 

where ¥ is the volume of propellant stored in the accumulator. 

The accumulator compliance a T is constant only for a linear accumulator; however, 
it is relatively simple to program a* as a function of line pressure for known nonlinear 
behavior such as adiabatic or isothermal compression of a gas. If a* is assumed to re- 
main constant and equal to its initial value over an incremental time period, small changes 
in volume of propellant stored in the accumulator can be related to small changes in pres- 
sure head as shown in the following equation: 

Ar = AH(a’) H==Hi {10) 

where Hj is the pressure head at the accumulator at the beginning of the interval. 

In order to properly simulate the cavitation -induced compliance at the pump inlet by 
the use of equation (10), the functional form of the pressure -volume relation for a volume 
of liquid containing vapor bubbles must be known. This relation, however, is not pre- 
cisely known. It would be expected to be nonlinear, and a hysteresis effect may be pres- 
ent because of a small nonuniform time lag related to heat transfer within the mixture 
during phase change. For some of the examples in this report the pres sure -volume curve 
for the fluid at the pump inlet will be approximated by assuming it to be hyperbolic in 
form and to have no time lag. This function is identical in form to the relation for iso- 
thermal gas compression in an accumulator. 

The compliance for an isothermal accumulator is given by 

- dr H o*-q 
dH h 2 

where H q and are the initial pressure and volume, respectively, of the gas in the 
accumulator. It should be pointed out that for a given compliance simulation (value of a') 
•f of the accumulator gas is considerably larger than the total volume of vapor cavities 
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Figure 5. - Conditions at pump before and after wave action. 


that would actually exist in the pump inlet region at pressure head H Q , because the phase 
change from vapor to liquid (or vice versa) associated with a change in head is an impor- 
tant contribution to the net incremental volume change of the fluid. 

Figure 5 shows conditions before and after wave action at a pump that is moving par- 
allel to the feed line with average velocity V p . Conditions are assumed to change at the 
instant the impinging waves leave the pump and then to remain unchanged for a time At. 
At that time new waves may impinge, and conditions may be once more changed. 

The continuity equation for the time period At expresses the condition that the flow 
past boundary AA (fig. 6) equals the storage in the accumulator plus the flow out of the 
pump (past boundary BB): 


< V 11 - V A f At - XT = V 22 A d At 

Introducing the relation for the accumulator (eq. (10)) and rearranging terms give 


v n - V P - 




in) 


The momentum equations are identical with the ones derived for the tank outlet (eqs. 
(2) and (3)). 

The head discharge -flow- rate relation for the pump after wave action is 

H 22 ' H ll = H 2 " H 1 + AH 2 + AH 22 ‘ AH 1 " AH ll = A + BQ + CQ 2 (12) 

The pump flow rate Q is measured with respect to the pump. 

Solving equations (2), (3), (10), (11), and (12) simultaneously for the discharge after 
wave action in terms of known conditions gives: 
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( 13 ) 


(■ 


2 f C 2 C l\ C i r 

C Q^ + B - — }Q + A + p - a - — - — = 0 


gA d gA t rj 


yg 


where 


a = H 2 + 2 AHg 


C 2 V 2 


0 = Hj + 2 AHj + 


c 1 v 1 


y = 1 + 


a'C 1 

gAtA f 


r = 


a' 

AtA, 


°iV|i 

,2iH i + ^ i ) tV p 


The positive root of equation (13) gives the flow rate through the pump after wave 
action. 

The resulting velocity in the feed line is given by 

v “-;(fr r ) (14) 

and the velocity in the discharge line is given by 

V 22 =-9- (15) 

A d 

Pressure wave magnitudes and line pressures are computed from equations (2)^ (3) ? 
(4), and (5). 

Injector . - The discharge line terminates with an injector, which forces propellant 
into the thrust chamber. 

The relation between flow through the injector and the net head across the injector is 
closely approximated by a square law relation (ref. 4) . This relation after wave action 
is 
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( 16 ) 


Q 2 2 

^=b i 2 ( h 11 -h 22 ) 

A d 

The chamber pressure is a function of flow through the injector. They are related by 

a combustion dead time r and a gas residence time 0_. The relation between chamber 

& 

pressure head H c and injector flow rate (ref. 1) can be written in the form 


dH H 
dt 


A T g6 g 


Q t (t - T> 


where c* is the characteristic velocity, A^, is the throat area of the engine nozzle, 
and Q.(t - r) is the injector flow rate at a time t-r in the past. In first order finite 
difference and wave-plan notation, equation (17) becomes 


H 22 " H 2 


+ xH 99 = yQ. 


where x = l/0 g , y ■ c*/A^,g0 g , At is the time increment between successive changes in 
the wave -plan analysis, and Q b is the burning rate of the propellant (Q b = Qj(t - t)). A 
higher order backward finite difference approximation may be employed with little addi- 
tional complication. 

The injector is located in a dome that has compliance. The net compliance of the 
dome region is somewhat larger than the discharge line compliance. This locally high 
compliance is also taken into account analytically through the use of a hypothetical accu- 
mulator. The mathematical expression used is identical with the equation for the pump 
inlet accumulator (eq. (10)). In a real system using hydrocarbon fuel much of the com- 
pliance in this region is related to structural springiness; therefore, a linear accumu- 
lator will be used in this area (a 1 = con- 
h a stant). For hydrogen -fueled engines with 

^ I V H 

v ! / v ll* n ll cavitation occurring on the dome side of 

iH i / AH U / ^H c 

the injector because of heating, a nonlin- 
v ’ ^ —PL/'/’S ear spring constant would be preferable. 

~~ ‘ )— / ) — / Figure 7 shows conditions before and 

*-*■ / \ - J j ( \ after wave action at the injector. The con- 

j tinuity equation expressing the condition 

a that the flow past boundary AA equals the 

Figure?. - Conditions at injector before end after wave action. flow through the injector plus the storage 


n 



for any small time At is 


Q 2 = V ll A d “ 


ATT 

At 


(19) 


Substituting the relation for the accumulator pressure -volume curve (eq. (10)) and de- 
noting the value of the slope as a** give 


«2 ‘ v n A d - < H n - H i>^' 


( 20 ) 


The line pressure after wave action is 


H 11 = H 1 + AH 1 + AH U = H i + 2 AH i + ~( y i " v n> 

g 


( 21 ) 


Combining equations (16), (18), (20), and (21) and solving for the line velocity after 
wave action give the following equation, which is quadratic in V^: 


“Ml + “ 2a V 11 + ^ ^ * H 22> = 0 


( 22 ) 


where 


*22 = 


yQ 1 + (H 2 /At) 
x + (1/At) 


Or = 1 + 


a"C 
g AtA. 


CV 1 

y s H, + 2 AH- + — 

g 


0 s L m 1 + 

A d At V 1 g / 


The positive root of equation (22) gives the discharge line velocity after wave action. 
The resulting injector flow rate is computed from equation (20). 
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Pump inlet perturbation suppressor . - A pump inlet perturbation suppression device 
is indicated in the model in figure 1 (p. 3) because of current interest in fix devices for 
reducing longitudinal structural oscillations of heavily loaded space launch vehicles. The 
function of such a device is to minimize the effect of flow perturbations induced in the feed 
system by structural motion. Two of the more common suppression devices being con- 
sidered are suction line volume compensating devices (ref. 6) and high compliance de- 
vices. 

The high compliance device can be represented by equation (10) by defining an appro- 
priate function for the parameter a' to fit the particular type of hardware (or pressure - 
volume relation) chosen. The volume flow into the device must be included as an addi- 
tional term in the continuity equation at the pump inlet (eq. (11)). For the special case 
where the losses and the inertia of the liquid in the standpipe are negligible and the com- 
pliance of the fix device is assumed to be of the same functional form as the cavitation in- 
duced pump inlet compliance, the net effect can be calculated with a single hypothetical 
accumulator by using the total compliance. 

The volume -compensating device (e. g . , a constant- volume bellows) can be repre- 
sented analytically by a similar approach by adding to the continuity equation the follow- 
ing term for the volume flow Q into the compensator (ref. 6): 

‘ 2S > 


where V p is the pump motion velocity and is the velocity of the other structural 
attachment point of the compensator. This point must be located in the forward part of 
the structure in order to get a large enough structural motion difference to operate the 
compensator. The kinematic factor k depends on the particular geometry of the com- 
pensator. Values of k for various types of compensators are given in reference 6. 
When the flow into the compensator is added to equation (11), the continuity equation be- 
comes 


V 
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H l ) = V 



(24) 


Equation (24) is used instead of equation (11) for determining the form of equation (13) 
when calculating the pump flow rate Q when a volume compensating fix device is present. 

Friction orifices . - Nonlinear (square law) orifices can be inserted throughout the 
system to simulate viscous effects. The analysis of this type of orifice subjected to pres- 
sure pulses is presented in reference 8. 
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Throttle valve . - A nonlinear (square law) throttle valve is located in the pump dis- 
charge line. The relation between pressure and flow is the same as that for a friction 
orifice and for an injector (eq. (16)). To throttle, the orifice coefficient is varied 
with time in a prescribed manner. For simplicity, the orifice coefficient is varied lin- 
early with time in the examples presented in this report. 

Digital Programming 

Time and position subscripts . - Combining the feed system component analyses into 
a single digital analysis of a propellant supply system is accomplished by introducing time 
and position subscripts. Through the use of these subscripts the time -independent equa- 
tions developed for a point are extended to time -dependent system equations. (A detailed 
explanation of the use of the subscripts is included in ref. 8. ) The subscripts are defined 
as follows: 

L position subscript 
J time subscript (time = J At) 

hi this notation the pressure head to the right of component L at time t - J At is writ- 
ten as HR(L, J). 

Each of the component analyses are written as a digital computer subroutine, which 
has the function of computing conditions at a system component at one instant of time in 
terms of conditions a small incremental time earlier, magnitudes of impinging waves, 
and occurrences that happend during the small time interval. The subscripted notation 
is generalized by defining variables as follows: 

HR pressure head to the right of the component 
HL pressure head to the left of the component 
AHR wave to the right and going away from the component 
AHL wave to the left and going away from the component 

If the component introduces no change of area, the liquid velocity in the line does not 
change and is defined as V. If the line area changes at the component, the velocity will 
change across the component and is defined as 
VI line velocity into component 
VO line velocity out of component 

The digital analysis consists of a finite difference type solution in which occurrences 
are computed at the end of small equal time intervals. The time interval is so chosen 
that wave travel times between adjacent components are always integer multiples of this 
time interval. The time subscript J refers to the integer number of time intervals 
since the initiation of the computations. The computations start with a known set of con- 
ditions at a discrete set of points and computes new conditions for these points for as long 
as desired. 
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The equations developed to analyze conditions at a component are directly applicable 
in the overall system analysis if the following substitutions are initially made: 

HL(L,J-1) 

h 2 = HR(L,J-1) 

AH^ = AHR{L-1, J-KX) 

AH 2 = AHLUh-1, J-KY) 

V 1= V<L,J-1) 

or 


Vj = VT(L, J-l) 

V 2 = VO(L,J-l) 

where the quantities KX and KY refer to the number of time intervals necessary for a 
sonic disturbance to travel from the respective adjacent components. 

The appropriate system component equations are then employed, and the corre- 
sponding outputs are as follows: 


H 11 = HL(L, J) 
H 22 = HR(L,J) 
AHjj = AHL(L, J) 
AH 22 = AHE(L, J) 
V n « VI(L,J) 
V 22 = VO(L,J) 


V n = V(L,J) 
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The preceding transformations result in a set of general relations that compute con- 
ditions at system component L at time J in terms of conditions that were computed in 
a previous cycle. 

Computer program s. - Digital computer subroutines were written for the principal 
system components. Subroutines were written for the tank outlet, the pump inlet station, 
the injector station, and the friction orifices. The throttle valve is a special case of the 
friction orifice subroutine in which the orifice coefficient is varied with time in a pre- 
scribed manner. A subroutine was not written for the fix device because in the examples 
of this report (see appendix B for system parameters) a high compliance device having 
the same functional form analytically as the pump inlet cavitation compliance was used. 
Thus the additional compliance could be incorporated into the accumulator at the pump in- 
let station. The four major subroutines are- given in appendix C along with examples of 
the FORTRAN programs. These subroutines are incorporated into the computer program 
for the monopropellant feed system. A copy of a typical FORTRAN IV program for the 
propellant feed system model is given in appendix D. 

Storage requirem ents. - Because of the large number of individual computations that 
must be made in the wave -plan analysis, it might be assumed that serious storage prob- 
lems would be encountered when computing solutions for complex propellant systems. 

This, however, is not the case. The digital computations associated with the wave -plan 
analysis may be carried out (starting with initial conditions) to any desired time. Be- 
cause the equations for each discontinuity are functions only of occurrences at the adja- 
cent discontinuities, only the information for the maximum number of time intervals re- 
quired for a sonic disturbance to travel from the farthest of the adjacent discontinuities 
need be stored. 

At the injector the additional requirement exists that information must be stored for 
at least the number of time intervals equivalent to the dead time. After these require- 
ments have been satisfied, the system conditions may be automatically reinitialized, and 
all necessary information will be available to the computer. Reinitializing allows a com- 
plex system composed of many discontinuities to be solved with relatively small computer 
storage capacity. 


RESULTS AND DISCUSSION 

The response of the hypothetical propellant feed system of figure 1 (p. 3) was compu- 
ted in order to illustrate the capabilities of the wave -plan method. Transient and steady - 
state solutions were determined for both periodic and ramp input disturbances. The nu- 
merical values of the feed system parameters used in the example are given in appen- 
dix B. The major parameters are summarized in table I. 
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Low-Frequency Periodic Inputs 

The response of the feed system was 
determined for low-frequency (1 to 16 
cps) sinusoidal inputs in the form of lon- 
gitudinal oscillatory structural motion at 
the pump and tank inlet. This response 
is representative of the type of low- 
frequency disturbances sometimes en- 
countered with spacecraft launch vehicles 
during the boost phase of flight. For 
simplicity the amplitude of the longitu- 
dinal structural oscillations was assumed 
to be independent of the magnitude of the 
thrust perturbations during the calcula- 
tion for each frequency. To simplify the 
example further, the amplitudes of the 
imposed structural velocity perturbations 
at the pump and tank were held constant 
over the frequency range studied. To 
establish a basis for comparison, the 
first series of calculations were made 
with no accumulators (zero additional 
compliance) at the pump inlet or at the 
injector dome. The first series of re- 
sults presented therefore pertain to a system in which the only compliance present is that 
due to the elasticity of the propellant and conduits. 

Frequency dependence . - The frequency response of the propellant system was de- 
termined. Calculations were made for a range of frequencies from 5 to 16 cycles per 
second. This range included the quarter-wave frequency for the 40-foot feed line, which 
for this case is 12. 5 cycles per second. The results are given in table U. All pressure 
heads were computed in feet of propellant. 

The system responses were generally periodic functions resembling sinusoids. Some 
slight nonlinearities were evident at lower frequencies. Figure 8 shows several typical 
periodic responses. * The magnitudes of the negative and positive perturbations are ap- 
proximately equal. 

*The flow and pressure head scales used in figures 8, 10, 12, 17, 18, and 20 are 
arbitrary with the zero corresponding to the local mean value or to the original steady- 
state value of the parameter. 


TABLE I. - MAJOR FEED SYSTEM PARAMETERS 


Line lengths, ft 


Feed line 

40 

Discharge line 

6 

2 

Line areas, ft 


Feed line 

0.2 

Discharge line 

0. 05 

Acoustic velocities, ft/sec 


Tank 

1000 

Feed line 

2000 

Discharge line 

3000 

Propellant density, lb/ft^ 

53 

Mean flow rate, ft /sec 

4 

Mean velocity in feed line, ft /sec 

20 

Engine parameters, msec 


Combustion dead time 

3 

Gas residence time 

i 

Pump coefficients 


A 

1208 

B 

668 

C 

-100 

Structural motion input velocity amplitude, ft/sec 


Pump 

0.8 

Tank 

0, 32 
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TABLE H. - AMPLITUDES OP PRESSURE AND FLOW PERTURBATIONS IN FEED SYSTEM 


(No accumulators in feed system (a' = 0, and a" = 0).] 


Frequency, 

cps 

Pressure head at feed 
line stations, a ft 

Pump inlet 
pressure 
head, 

- ft 

Pump outlet 
pressure 
head, 
ft 

Injector 
flow rate, 
ft 3 /sec 

Qi 

V BB A f 

Dome 

pressure 

head, 

ft 

Chamber 

pressure 

head, 

ft 

2 

3 

4 

5 

16,1 

23. 1 

28,8 

35, 9 

28.4 

0,037 

0.23 i 

25.7 

10.9 

7 

26.2 

37.3 

47,4 

55.9 

44.2 

,057 

.36 

39.2 

16. 9 

9 

41. 9 

58.3 

73.4 

82. 7 

65.5 

.084 

,52 

59.8 

25.2 

10 

52.7 

73, 6 

90,1 

101, 2 

79,5 

, 103 

,64 

72, 9 

30,8 

11 

66.7 

93,3 

111.9 

123,1 

96, 1 

, 126 

.79 

88.2 

37,3 

12 

80.8 

112. 2 

132,7 

142. 0 

110, 7 

,145 

,91 

101.8 

42.9 

12, 5 

87.5 

120. 5 

141.4 

148.5 

115. 5 

,151 

.94 

106.4 : 

44,9 

13 

93.6 

123.8 

148.0 

149,7 

117.6 

. 152 

.95 

107.9 

45,4 

13.5 

94.1 

126.4 

146.9 

148,5 

116.7 

.150 

* 94 

106. 1 

44,7 

14 

93.1 

124.7 

142.1 

142.2 

111.2 

,144 

* 90 

101.5 

42.8 

14,5 1 

88.8 ; 

122.2 

134.2 

133,8 

105,0 

,135 

,84 

94.6 

40,2 

16 

76.9 

101.9 

108. 7 

99,7 

78,9 

.102 

.64 

70,7 

30, 1 


a Feed line stations 2, 3, and 4 are located 24, 16, and 8 ft above the pump inlet station. 


Pump suction pressure . - The performance of the pump is directly related to the in- 
let (suction) pressure. Because of this relation, pressure fluctuations at the pump are 
of prime importance, and the pump must be designed to operate satisfactorily under the 
most extreme pressure perturbation anticipated. As can be seen in table II, the maxi- 
mum pressure perturbation (149. 7 ft, 55 psi) occurred slightly above 13 cycles per sec- 
ond, which is near, but somewhat greater than the quarter -wave frequency (12. 5 cps) for 
the 40-foot feed line. The frequency dependence of pump suction pressure perturbation 
amplitude can also be obtained from table II. 

Minimum feed line pressure . - Relatively low pressures exist in the pump suction 
line. Therefore the lowest pressure that can occur in the feed line must be determined 
if the occurrence of feed- line cavitation is to be predicted. The instantaneous pressure 
at any point in the feed line is given by the sum of the steady-state or mean value and the 
perturbation value at that instant. The location in the feed line where this sum is a mini- 
mum must be determined, and the magnitude of the minimum value must be calculated. 
For the frequency range studied (table n), the maximum pressure perturbation amplitude 
occurs at the pump inlet, and the size of the perturbation amplitude decreases with dis- 
tance up the feed line. 

At any instant the local steady-state (or mean) pressure occurring in a vertical pipe 
line in a uniform acceleration field increases linearly from the tank to the pump. The 
pressure head gradient (feet of head per linear foot) due to acceleration of the center of 
mass of the vehicle in a gravitational field is 
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Figure 9. - instantaneous minimum values of local pressure in feed 
line. 


dHg a z + g 
dz g 

where a 2 is the vertical acceleration of the center of mass of the vehicle. The pressure 
gradient due to friction losses is given by the Darcy equation as 

dHj = _„2 

dz 2gD 

Thus the overall steady (or mean) pressure gradient is 

dH _ a z * * fV 2 
dz g 2gD 

For the system studied the pressure gradient due to friction is relatively small compared 
with the gravitational gradient and has a value of 0. 4 foot of head per linear foot. 

Figure 9 shows the variation in the local instantaneous minimum pressure along the 
feed line computed for two arbitrarily chosen overall steady gradients. The variation is 
given for three different structural oscillation input frequencies. The dashed line is for 
a slope of 2.0, which corresponds to an acceleration of the center of mass of 1. 4 g's. 

The solid line is for a pressure slope erf 3. 0 corresponding to an acceleration of 2. 4 g's. 
The lowest pressure generally occurs between 8 and 24 feet from the pump inlet and is 
considerably lower than the minimum value of pressure at the pump inlet. 
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Feed-tine flow rat?, cu ftlsec Pump flaw rate, cu ftisec Chamber pressure head, ft Pump suction pressure head, ft 




This analysis points out that the location in the feed line of maximum negative pres- 
sure perturbation is not necessarily the point experiencing the lowest instantaneous net 
pressure. 

Injector flow rate perturbation. - The injector flow rate perturbation is of particular 
importance because longitudinal structural oscillations of a launch vehicle are driven by 
thrust perturbations, which in turn are directly related to the propellant flow rate through 
the injector into the thrust chamber. As shown in table II, the frequency dependence of 
the amplitude of the propellant flow rate perturbation through the injector follows a pat- 
tern similar to the pump inlet pressure perturbation. 

The amplitude of the ratio of the injector flow rate perturbation to the equivalent flow 
rate perturbation caused by the pump motion is also given in table II. This ratio reaches 
a maximum of 0. 95 at the equivalent quarter- wave resonance condition (13 cps) of the feed 
system. At quarter-wave resonance the impedance to perturbation flow up the feed line 
toward the tank is very high. This high value indicates that as the pump moves up and 
down only a small part (5 percent) of the flow disturbance generated can travel up the line 
to the tank. Instead most (95 percent) of the disturbance travels through the pump and in- 
to the engine where it results in thrust perturbations. 


Effect of Pump and Injector Compliance on Low-Frequency Oscillations 

The accumulators (fig. 1, p. 3) at the pump inlet and injector were used to simulate 
locally high compliances at these stations. The effect of these compliances on pump inlet 
pressure perturbations, pump inlet flow perturbations (measured in the feed line just 
above the accumulator simulating inlet compliance), pump outlet flow perturbations, and 
thrust chamber pressure perturbations were calculated with the model. 

Pump inlet with linear accumulator. - A quantitative indication of the effect of 

cavitation-induced compliance at the pump inlet was obtained with a constant compliance 

-5 -5 

accumulator. The values used in the calculations were a' - 0. 2x10 , 1. 0x10 , and 

2. 0x10 cubic feet per foot head. 

Figure 10 shows the computed responses of various system parameters to sinusoidal 
longitudinal structural oscillations (compare with fig. 8). The calculations were carried 
out for a frequency of 10 cycles per second and an inlet compliance of 2x10”^ cubic feet 
per foot. The resulting perturbations in pump suction pressure, thrust chamber pressure, 
pump flow rate, and feed line flow rate are shown. The most noticeable effect of the ac- 
cumulator is that system nonlinearities originating in the tank become apparent in the 
feed line flow rate. The nonlinear disturbance is transient in nature and eventually dies 
out. The frequency dependence of the amplitude of the pump inlet pressure perturbations 
is shown in figure 11 for the three compliance constants. The response curve with no ac- 
cumulator in the line is also given for comparison. The general trend is for the resonant 
22 
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Figure 11. - Frequency response of pump suction pressure pertur- 
bation with iinear accumulator at pump iniet. 

frequency to decrease as the inlet compliance is increased. Thus cavitation induced com- 
pliance in the pump inlet region tends to lower the resonant frequency of the suction line. 
Other parameters such as injector flow rate perturbation and thrust chamber pressure 
perturbation are affected in a similar manner. 

Pump i nlet with hype rbo lic accumulator . - An indication of the effect of nonlinear 
compliance at the pump inlet was obtained by calculating the response of the system with 
the hyperbolic accumulator. This type of accumulator may be used to approximate pump 
inlet cavitation, a M high compliance type" fix device, or a combination of both. 

Accumulators having the characteristics shown in the table above were analyzed. 

The smaller initial volumes can be considered to be representative of the cavitation- 
induced inlet compliance at pump suction pressures greater than or equal to H C1 (fig, 4, 
p. 7), whereas the larger initial volumes (0. 01 and 0, 1 ft ) are representative of gas ac- 
cumulator fix devices or alternately a combination of both compliances. The larger ini- 
tial volumes are also typical of cavitation compliance for pumps operating at suction 
pressures lower than Hgj (fig. 4). 

The nonlinear response of the system with hyperbolic accumulators at the pump inlet 
is indicated by the pressure -time and flow- rate -time traces shown in figure 12. The re- 
sponses for system parameters at two arbitrary sets of conditions are shown. For the 
majority of cases calculated the responses of the pump suction head, chamber head, and 
pump flow rate are characterized by sharp high peaks and flat valleys. Occasionally the 
traces were quite irregular because of the presence of higher harmonics at a particular 
condition (near resonance). It was also noted that very large accumulator compliances 
tended to lengthen the time for the starting transient to decay. The feed line flow rate 
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Frequency, cps 


Figure 13, -Nonlinear response of pump suction pressure with gas 
accumulator. 


above the accumulator displayed no 
predominant pattern but was charac- 
terized by nonlinear behavior. 

The frequency response of pump 
inlet pressure perturbation amplitude 
for a pump having a hyperbolic accu- 
mulator at the inlet is compared in 
figure 13 with the response of a pump 
with a linear (or constant compliance) 
accumulator at the inlet. The hyper- 
bolic accumulator has the same com- 
pliance (2xl0~^ ft 3 /ft) at mean line 
conditions as the linear accumulator. 
The hyper l)o lie accumulator is repre- 
sented by an isothermal gas accumu- 
lator with a volume of 0. 004 cubic feet 
(6. 9 in. ) at a mean line pressure of 
200 feet. The curves for the positive 
and the negative amplitudes as a func- 
tion of frequency are mirror images 
for the linear accumulator. The curves 


for the hyperbolic accumulator, however, are not mirror images because, as shown in 
figure 12, the positive amplitude is greater than the negative amplitude. This results 
from the hyperbolic form of the compliance, which increases with decreasing pressure 
and decreases with increasing pressure. 

The effect of pump inlet compliance on the suction line resonant frequency was de- 
termined by examining the computed system response with hyperbolic accumulators of 
various sizes at the pump inlet. In figure 14, the positive and negative amplitudes of the 
perturbation pressures are plotted as functions of frequency for accumulators of various 
sizes. The accumulators have the various indicated volumes at an initial pressure of 
150 feet. Figure 14 shows that increased pump inlet compliance has the effect erf lower- 
ing the suction line resonant frequencies, as was noted for the linear accumulator results 
(fig. 11, p. 23). Because of the hyperbolic functional form of the accumulators, the posi- 
tive amplitudes at resonance are greater than the resonant value for the system with no 
accumulator. Only at very large compliances does the positive amplitude at resonance 
become equal to or smaller than the resonant value for the system with no accumulator. 
The negative amplitude at resonance is always smaller than the no-compliance case and 
decreases steadily with increasing values of inlet compliance. The peak to peak ampli- 
tude at resonance first increases with the addition of compliance (V Q = 0. 0005) then 
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Frequency, cps 
\b) Negative amplitudes. 

Figure 14, - Pump suction pressure parturition with hyperbolic accumulator. 


steadily decreases with increasing compliance. Figure 14 demonstrates that any increase 
in pump cavitation occurring during the boost phase of flight would lower the resonant 
frequency of the feed system. Further, a pure compliance (no inertia, or resistance) 
fix device placed at the pump inlet would have the function of lowering the pump suction 
line resonant frequency. 

Injector dome co mp li ance . - The propellant feed system response was determined 
for three values of compliance of the dome. The compliances are assumed to remain 
constant, and the values used include 5X10 2. 5x10 and 0. 5x10 ® cubic feet per foot. 
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Figure 15. - Frequency dependence of pump suction pressure per- 
turtation with accumulator at injector dome. 


The frequency dependence of the pump inlet pressure perturbation is shown in fig- 
ure 15. The curve for a compliance of 0. 5x10 is not shown because it follows the zero 
compliance curve so closely. This value is close to the value for the line compliance. 

The accumulator at this point in the system has no great effect at longitudinal structural 
oscillation frequencies. The critical frequency and maximum amplitudes are both slightly 
lowered. The accumulator has a relatively small effect at this location because the lin- 
earized value of the injector impedance (at the mean flow condition) is only about one-half 
of the characteristic impedance of the discharge line. 


Transient Response of System 


Throttling . - The digital model was employed to determine the propellant system re- 
sponse to throttling. Throttling transients were introduced at the pump discharge line 



Figure 1 a, * Variation of orifice coefficient with time during 
throttle operation. 


throttle valve (fig. 1, p. 3). With the sys- 
tem running at a steady- state operating con- 
dition the throttling transient was introduced 
by initiating a linear variation of the orifice 
coefficient from B^ to B t2 in a short 
period. This variation is illustrated in fig- 
ure 16. 

When the system was subjected to tran- 
sient inputs due to throttling it responded 
in one of two ways. The system flow and 
pressure usually underwent short duration 
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disturbances that eventually died out. The system then returned to the original steady- 
state values or to a new steady-state condition when a permanent change in the system 
operating conditions had been introduced. Alternately there were conditions under which 
the initial throttling transient disturbance did not fully damp out before the system went 
into an unstable mode of operation (chugging) because it ended up in an unstable operating 
region. 

Figures 17 and 18 contain examples of throttling transients. In figure 17 the system 
was throttled while oscillating longitudinally. (No lateral motions were assumed to be 
present. ) Both the pump suction pressure and the thrust chamber pressure perturbations 
were approximately steady- state sinusoids before the throttling was initiated. Following 
the short transient period, the pump suction and thrust chamber pressures quickly re- 
turned to a stable steady -state sinusoidal oscillation at a new mean level and with a new 
perturbation amplitude. The amplitude of the structural motion input was held constant 
during the computational period. The change in perturbation amplitudes was caused by a 
shift in system operating conditions. 

Figure 18 shows responses for a system that was throttled from an initially steady 
operating condition (no structural oscillations). In figure 18(a) the throttling introduced a 
transient that quickly died out, and the system ended at a new stable operating point. Fig- 
ure 18(b) shows a condition where the system was throttled to a greater degree. In this 
case the system mean operating parameters changed significantly, and the system went 
into an unstable high frequency chugging oscillation before the throttling transient had 
died out. The calculations were terminated at the end of the pressure -time traces shown 
because the system was unstable. 

Girobal snubbing. - During flight one or more engines of a rocket may use gimbaling 
action to keep the missile on the desired flight path. When the thrust chamber and feed 
system pumps are mounted so that they both move during the gimbaling action, transients 
are introduced into the system. During the initial gimbaling maneuver, the accelerations 
are low because of engine inertia, and no noticeable transients would be introduced into 
the feed system; however, under some circumstances the engine gimbal motion may be 
stopped suddenly by snubbing action. A sudden stop would introduce a transient into the 

feed system. In this example the pump-engine assembly 
will initially be assumed to be moving relative to the tank 
outlet station at the maximum velocity (2 ft /sec) obtained 
duri ng the gimbaling motion. The motion of the engine - 
pump assembly is then assumed to be linearly snubbed 
from 2 feet per second to rest in 40 milliseconds (fig. 19). 
The effect of the snubbing action on the pump suction pres- 
sure head is shown in figure 20. The snubbing action was 
initiated at 10 milliseconds, and the system response was 



Figure 19. - Pump velocity during 
gimbaling ’Snubbing maneuver. 
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Figure 20, - System response to rapid snubbing of pump motiw resulting from engine gimbaling. 


calculated out to 600 milliseconds. Figure 20 shows that in this severe case snubbing 
caused an increase in pressure of 97 feet above and 41 feet below the mean conch tion. 
The pressure transient quickly dies out, however. 


CONCLUDING REMARKS 

The simple propellant feed system analyzed in this report does not correspond to any 
operational propellant supply system. It is, however, fundamentally similar in many 
respects to operational systems and serves to illustrate how the wave- plan method can be 
applied to analyze unsteady flow in a rocket propellant feed system. 

The fact that a monopropellant (or one side of a bipropellant) feed system was ana- 
lyzed is not meant to imply that the method is limited to such a system. The method can 
readily be extended to a bipropellant system by coupling two feed systems together at the 
thrust chamber. This coupling is done through the relation between thrust chamber pres- 
sure and the fuel and oxidizer weight flow rates; however, for some analytic purposes the 
cross -coupling effects on the individual feed system responses are small, and the simpler 
monopropellant model can be used to explore the effect of varying basic system param- 
eters and geometries. 

In the calculated examples for the feed system response to structurally imposed per- 
turbations, the amplitude of the structural motion at any input point was assumed to re- 
main constant during the calculation; that is, no feedback was provided between combus- 
tion chamber pressure (and, therefore, thrust) perturbations and the structural oscilla- 
tions. Similarly no feedback was provided between local feed system pressures (at the 
pump inlet for example) and the structure. However, the model can be readily combined 
with structural dynamic equations to provide the necessary feedback. Thus, in cases 
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where the structural dynamics equations are available the wave -plan method of analyzing 
the feed system can be incorporated into a nonlinear stability analysis. 

Because of the present lack of dynamic data on turbopumps the steady- state operating 
curve was used to determine the dynamic resistance of the pump. In addition, the analysis 
was limited to the design operating region of the pump. Again, these assumptions were 
not basic to the method. A very complex pressure -flow map or tabulated dynamic data 
for the pump can be easily read into the digital program when available. Further, no 
coupling between the pump and turbine was put into the examples. Including the turbine 
would require a method of characteristics solution on the gas side as well. This exten- 
sion, although entirely feasible, was considered beyond the scope of this report. 

In conclusion it should be emphasized that the calculated response of a fluid system 
is very specific and highly dependent on the particular hydraulic and structural dynamic 
characteristics of the system analyzed. It is, therefore, very difficult and usually 
dangerous to make conclusions that apply to systems in general. The best approach is to 
analyze the specific hardware at the operating conditions of interest. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 23, 1966. 
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APPENDIX A 


SYMBOLS 


A, B, C 

coefficients for static net head- 
discharge curve for pump 

C 2 

wave velocity in line following 
component, ft/sec 

a, b, c 

coefficients for quadratic 

C* 

characteristic velocity, ft /sec 


functions 

D 

line diameter, ft 

A d 

cross-sectional area of dis- 
charge line, ft 2 

f 

g 

friction factor 
acceleration due to gravity, 

A f 

Arp 

cross-sectional area of feed 
line, ft 2 

throat area of engine nozzle. 

H 

ft/sec 2 

pressure head, ft 

X 

9 

fr 

cross-sectional area of tank. 

AH 

pressure head change across 
wave or component, ft 

ft 2 

H c 

chamber pressure head, ft 

a' 

compliance of pump inlet 

Q 

accumulator, ft /ft 

H f 

pressure head due to friction 
losses, ft 

a M 

compliance of injector dome 
accumulator, ft /ft 

H g 

pressure head due to accelera- 
tion of gravity, ft 

b f 

frictional orifice coefficient, 
2 1/2 
(ft/seO 

H o 

H 1 

initial pressure head at 
accumulator, ft 

pressure head to left of com- 

B I 

orifice coefficient for injector, 
9 1/2 
(ft /sec 4 ) 


ponent before wave action, 
ft 

b t 

C d 

orifice coefficient for throttle 
2 ^/ 2 

valve, (ft /sec ) 
wave velocity in discharge 

AHj 

wave impinging from left of 
component before wave 
action, ft 

Q 

line, ft/sec 

H 2 

pressure head to light of com- 
ponent before wave action. 

c f 

wave velocity in feed line, 
ft/sec 

ak 2 

ft 

wave impinging from right of 

c t 

wave velocity in tank, ft/sec 


component before wave 

c i 

wave velocity in line preceding 
component, ft/sec 


action, ft 
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H u 

pressure head to left of com- 
ponent after wave action, ft 

V BB 

amplitude of pump sinusoidal 
structural motion velocity, 

AH U 

wave leaving component on left 


ft/sec 

side after wave action, ft 

V i 

linear velocity of structure at 

h 22 

pressure head to right of com- 
ponent after wave action, ft 


compensator attachment 
point, ft/sec 

AH 22 

wave leaving component on 
right side after wave action, 

V 

P 

linear velocity of pump due to 
structural motion, ft/sec 


ft 

V t 

linear velocity of tank due to 

Hg, Hgj 

pump suction pressure head, 


structural motion, ft/sec 

AH P 

ft (see fig. 4) 

pressure head rise across 
pump, ft 

V 1 

velocity in line to left of com- 
ponent before wave action, 
ft/sec 

k 

compensator kinematic factor 

V 2 

velocity in line to right of com 
ponent before wave action, 

L d 

length of discharge line, ft 


it /sec 

L f 

length of feed line, ft 

V 11 

velocity in line to left of com- 

Q 

volume flow rate, ft^/sec 

nent after wave action. 


burning rate of propellant, 


ft/sec 

Q b = Qj(t - t), ft 3 /sec 

V 22 

velocity in line to right of com 

Q„ 

volume flow into compensator. 


ponent after wave action, 

c 

ft 3 /sec 


ft/sec 

Qj f Q| 

injector discharge at time 

Y 

volume, ft 


equal to dead time in past, 
ft 3 /sec 

r 0 

initial volume of gas in 

3 

accumulator, ft 

q 2 

injector discharge after wave 
action 

x,y 

injector delay and dead time 
coefficients 

R ! 

ratio of feed line to tank area 


gas residence time, sec 

t 

V 

time, sec 
velocity, ft/sec 

7 

dead time, sec 
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APPENDIX B 

SYSTEM PARAMETERS USED IN EXAMPLE 


The hypothetical propellant feed system is shown in figure 1 and consists of a tank 
and a pump connected by a feed line containing a non compensating bellows. The pump 
discharges into a discharge line and into a nonlinear (square law) injector. 

In order to make a digital analysis of the response of the system it is necessary to 
assign numerical values to all parameters. The effect of certain key parameters on the 
response of the system can be evaluated by varying their numerical values. 

The propellant tank and line parameters were assigned the following values: 


Wave velocity in tank, C t , ft/sec 

Wave velocity in feed line, C f , ft/sec . . . 

Wave velocity in discharge line, C^, ft/sec 

Length feed line, L f , ft 

Length discharge line, L d , ft 

Area tank, A*., ft^ 

£ 2 

Area feed line, Ap ft 

Area discharge line, A^, ft^ 


1000 
2000 
3000 
40 
. 6 
20 
0.2 
0.05 


The propellant was assumed to have a density of 53 pounds per cubic foot. The mean 
propellant flow rate was 4 cubic feet per second at the design operating condition. The 
injector pressure drop was approximately 40 percent of the dome pressure. An injector 
coefficient Bj equal to 2. 8 was used. 

Thrust chamber coupling with the feed system was provided by using a chamber pres- 
sure that is a function of the injector flow rate and related in time by a time delay factor 
(gas residence time) and a dead time. The combustion dead time was chosen to be 3 mil- 
liseconds. The gas residence time was taken as 1 millisecond, and the quantity 
c*/Aj.gf?g was approximately 370 000 per square foot. 

The following pump coefficients were employed in equation (9): A = 1208, B = 668, 
and C = -109. 

Viscous losses in the tank were very small and were neglected. The total friction 
loss in the feed line at mean line conditions (Q * 4 ft 3 /sec) was assumed to be 16 feet. 

The corresponding loss in the discharge line was assumed to be 100 feet. The higher 
value is due to the high velocities present in the discharge line. These losses were dis- 
tributed over four friction orifices located in the feed line and one located in the discharge 
line. 

At the pump and tank, coupling to the structure was provided by forced motion of 
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these components along the vehicle vertical axis. This structural motion was assumed to 
be sinusoidal, and the amplitude of the tank motion was 40 percent of the pump motion 
amplitude. The amplitude of the structural velocity perturbation at the pump was taken 
to be 0. 8 feet per second in the subsequent examples. This value is 4 percent of the mean 
flow velocity in the feed line. For simplicity, the same structural motion amplitudes 
were used at all frequencies. 

To complete the definition of mean operating conditions for the system the height of 
fluid in the tank and the pressure at some point in the system must be designated. Be- 
cause of the long propellant line length, the dynamic response of this system has only a 
small dependence on propellant height in the tank. Since critical frequencies usually 
occur near burnout, a constant level of fluid in the tank of 2 feet was assumed and was 
used for the subsequent calculations. The propellant level could have been treated as a 
time-dependent variable, but the calculations were carried out for only a short interval, 
and the added complexity was not warranted. The dynamic response of the system as 
formulated is independent of mean pressures except for cases when local compliances 
are pressure dependent. For these cases the mean pressure at the local station is de- 
fined in the RESULTS AND DISCUSSIONS section of the report. 
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APPENDIX C 


COMPUTER SUBROUTINES 

With the use of the subscripting notation four digital computer subroutines were 
written. These subroutines perform the following operations: 

(1) Equate the nonsubscripted initial conditions to their subscripted counterparts 

(2) Solve the appropriate equations for a new set of conditions 

(3) Equate the newly computed conditions to their subscripted counterparts 
The four principal subroutines are as follows: 

(1) Subroutine tank 

(2) Subroutine pump 

(3) Subroutine injector 

(4) Subroutine friction orifice 

Each subroutine requires the definition of the time counter J and the position 
counter L before it is entered. In addition, the following quantities must be defined for 
each subroutine: 

Subroutine Tank: 

VT average velocity of tank over time interval 

KX number of time increments to nearest component in tank 

KY number of time increments to nearest component in feed line 

Cl wave velocity in tank 

C2 wave velocity in feed line 

R ratio of feed line to tank area 

Subroutine Pump: 

VP average velocity of pump over time increment 

KX number of time intervals to nearest component in feed line 

KY number of time intervals to nearest component in discharge line 

Cl wave velocity in feed line 

C2 wave velocity in discharge line 

A, B, C pump characteristic coefficients 
AP compliance of pump inlet spring-piston accumulator 

R ratio of feed line to discharge line area 

Subroutine Injector: 

BI injector coefficient 

KY number of time intervals to nearest component in discharge line 
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C wave velocity in discharge line 

APP compliance of injector spring-piston accumulator 

Subroutine Friction Orifice: 

BF friction orifice coefficient 

C wave velocity in line 

K number of time intervals to adjacent component 

The FORTRAN IV program for each subroutine including tables relating subscripted to 
nonsubscripted variables is as follows: 
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SUBROUTINE TANK (V1»V2»H1 ,H2,DH1 »DH2 , VI 1 »V22 ,H1 1 »H22 ,DHII » DH22 »VT ) 
COMMON/B0X1/AP,APP»AD»AF,DT,X,Y»AT .80 
C0MM0N/B0X3/A,B,C»G»BF, ERROR 
COMMON /B0X2/C1 »C2 .R 


VI 

VI f L * J-l ) 

V2 


HI 

HL(L tJ-1 )' 

H2 

HR (L* J-l 3 

0H1 

W8 f L-l iJ-KX3 

DH 2 

OHLCL+lf J-KY) 

Vll 

VI ( L * J) 

V22 

V0IL*J) 

Hit 

HL ( L * J ) 

H22 

HRILfJy 

0H11 

DHL ( L # J ) 


DH22 OHR(L.J) 

VT VT<J> 

AA=(1.-R*R)/(2.*G> 

BB=C1*R/G+C2/G-VT*<R-R*R)/G 

CC=H2+2 •♦0H2-H1-2 •*0H1-C2*V2/G-C1*V1 /G+C1*VT* (I .-R ) /G-VT*VT* 
1 ( 1 • -R I * ( 1 ,-R ) / ( 2 • *G ) 

V22 * V2 
DO 511=1,20 

2 ER = (AA*V22*V22+BB*V22+CCl/t2.*AA*VP2+BB) 

V22 * V22-ER 
I F { ABS ( ER ) -ERROR ) 1.1,5 
5 CONTINUE 
£RROR= 1000* 

WRITEI6.100) V22.ER 
100 FORMAT ( 2E15.8 ) 

1 Vll * VT#(1.-R)+V22*R 
0H11 = DH1+C1/G*(V1-V11 ) 

DH22 * DH2+C2/G* I V22-V2 ) 

HI 1 = Hl+DHl +DH11 
H22 = H2+DH2+DH22 
RETURN 
END 


SUBROUTINE I N JECT (H1.H2.V1.01 ,DH1 ,H1 1 ,H22 ,V1 1 ,02 »DH11) 
COMMON /B0X1/AP,APP,AD,AF,DT ,X ,Y»AT ,BI 
COMMON /B0X3 / A.B.C.G.BF.ERROR 

DOUBLE PRECISION APPD .CD ,GD ,DTD . DH1D . V1D .HID , VI 1 D ,B I D.H22D . 
1 ALPHA, BETA, GAM, AA.BB.CC.ERD 
H2 HR(L.J-l) 

VI V (L.J-1) 

HI HLIL.J-1) 

01* 0(L,J-K5l 

DH1* DHR f L-l » 1 

VII V IL.JI 

HU HL(L.J) 

H22 MR ( L » J ) 

0H1 1 DHL(L.J) 

Q2 = Q ( L * J ) 

H22*(Y*01+H2/DT)/(X+1./DTJ 

APPD=A»P 

CD=C 

GD=G 

DTD=DT 

DH1 D=DHl 

V1D-V1 


41 



H10=H1 
B T D=B ! 

H22D=H22 

ALPHA=1.D0+APPD*CD/ (GD*DTD) 

BETA=APPD/DTD»(2.D0*DH1D+CD*V1D/GD) 

GAM«H1D+2.D0*DH1D+CD/GD*V1D 

VlID=VlD 

AA * ALPBA*ALPHA 

BB=BID*BID*CD/GD-2,00*AL B HA*BFTA 
CO-BIO*BID*(GAM-H22D)+BETA *BETA 
NN=0 

1 ERD=( AA*V11D*V11D+8B*V11D+CC ) / ( 2 . D0*AA*V1 1D+8B ) 
V11D=V1 10— ERD 

ER=ERD 

Vll=VllD 

IF (ABS(ER) -ERROR) 2*2*6 

6 NN=NN+1 

IF<NN~50) 1,1,7 

7 WRITE(6*101)V11 , ER 
ERROR = lOOO • 

101 FORMAT ( 2E16, S ) 

2 0H11 = 0H1+C/G*( Vl-Vll ) 

HI 1 = Hl+DHl+DHll 
Q2=V11*AD-APP*AD/DT*(H11-H1 ) 

RETURN 

END 


SUBROUTINE FOR { HI ,DH1 ,DH2 ,H2 , VI ,H1 1 ,DH22 , DH1 1 ,H22,V11 ) 
C SUBROUTINE FRICTION ORFICE 


c 

HI 

HL I L , J— 1 ) 

c 

DH 1 

DHRtL-1 

c 

D H2 

DHL (L+l , J-K ) 

c 

H2 

HR ( L » J~ 1 ) 

c 

VI 

V f L , J— 1 ) 

c 

HI 1 

HLCL, J) 

c 

DH 11 

DHLIL »J) 

c 

DN22 

DHRCL»J) 

c 

N22 

HRILjJS 

c 

Vll 

V ( L , J ) 


COMMON /B0X3 / A, B, C,G,BF, ERROR 
B8 = 2,*BF*BF*C/G 

CC=-(H1+2.*DH1+2.*C/G#V1-(H2+2.*DH2 ) )*BF*BF 
IF (CC) 3,3,2 

2 BB = -BB 
CC = -CC 

3 Vll = Vl 

DO 611=1,20 

A ER = <Vll*Vll+BB*VIl+CC)/(BB+2.*Vll> 

Vll = Vll -ER 

IF (ABS(ER)-ERROR ) 5,5,6 
6 CONTINUE 
ERROR* 1000 . 

WRITE <6, 100) V22 , ER 
100 FORMAT ( 2E1 5 • fi ) 

5 DH1 1 = DH1+C/G* ( Vl-Vl 1 ) 

HI 1 = Hl+DHI+DHll 
DH22 = DH2+C/G*<V11-V1 ) 

H22 = H2+DH2+0H22 

RETURN 

END 
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SUBROUTINE PUMP (VP ,V1 »V2 *H1 , H2 » DH1 ,DH2 »V11 ,V22 ,H11 »H22 »DH1 1 .DH22 ) 
DOUBLE PRECISION QD * AA ,BB ,CC .C2D ,BD,GO »C1D , AFD«APD»DTD,H2D» 
1DH2D,DH1D,V1D,V2D,AAD»VPD,ADD,ERD 
COMMON /B0X1/AP,APP ,AD*AF ,D? ,X ,Y , AT t BO 
COMMON /BOX2 /Cl »C2 »R 
COMMON /B0X3/ A »B.C,G,BF, ERROR 

VI * VI ( L » J— 1 ) 

V2 * V0(1,J-1) 

HI = HLIL,J-1) 

H2 = 

DH1 - DHR(L-1»J-KX) 

OH 2 = DHL ( L+l » J— KY ) 

VII a VI(L*J) 

V22 a VO(l.J) 

HI 1 = HUL.J) 

H22 * HR ( L » J I 

DH11 » DHL ( L » J ) 

DH22 = DHR(LiJ) 

0 = V2*AD 
H1D=H1 

C2D = C2 
BD s B 
GD = G 
CIO = Cl 
AFO = AF 
APD » AP 
OTO * DT 
H2D = H2 
DH20 = OH 2 
OHIO = DH1 
VI D = Vl 
V20 = V2 
AAD a A 
VPD = VP 
OD = O 
ADO = AD 
H22D = H22 
AA=C 

BB = BD-C2D/(GD*ADD)-ClD/(GD*AFO*f 1.DO+APD*C1D/(GD»DTO) n 
CC = -(H2D+2.D0#DH2D-H1D-2.D0*DH1D— (C1D*V1D+C2D*V2D» 
1/6D-AAD+C1D/GD*( APD/DTD* < 2 . DO*DHl D+C1D*V1D/GD J+VPD ) / 

2( 1.D0+APD*C1D/<GD*DTD> ) > 

DO 511*1,20 

2 ERD = ( AA*QD*GD + BB#QD+CC > / t 2 » D0+AA*GD+BB ) 

QD a QD-FRD 

ER = FRD 
Q = QD 

IF ( AB5(ER)-£RR0R ) 3,5,5 
5 CONTINUE 
£RR0R=1000. 

WRITE (6, 100) V22.ER 
100 FORMAT ( 2E 1 5 . 8 ) 

3 V22 = Q/AD 

Vll = ( (APD*(2.DO*DH1D+C1D*V1D/GDH/DTD+VPD+QO/AFD)/ 

1 (1.D0+APD*C1D/ (GD*DTD) ) 

DH1 l*DHl+Cl /G* t Vl—Vl 1 1 
DH22=DH2+C2/G*(V22~V2 ) 

HI 1 =H1 +DH1+DH1 1 

H22=H2+DH2+DH22 

RETURN 

END 
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APPENDIX D 


EXAMPLE OF PROGRAM FOR MONO PROPELLANT FEED SYSTEM 



Read 


Compute and assign initial conditions 

Start-* — 

program 


in tank, feed line, discharge tine and 


data 


at tank outlet, pump, and injector 



Free surface reflection 
DPR(t ( JI*-DPU2 ( J-K1) 



Tank friction 
through B 
L-£U1 ‘I 



R ^Sub routm*P \ 8 
frictlcm orifice/^" 




Feed line friction 
through C 
L p LI + 1, J p L2 - 1 

IT 

Discharge line fric- 
tion through D 
L-L2*l> L3-1 


VP - VBS sin 2nfJAt 
1-12 C2-CD 
R*R2 O-KI 
Cl-CF KY * &2 



10 s — ^ 

( Subroutine \ 
pump J~ 


B — — " — 14 

/ Subroutine \A_ 
l injector ) f 


i — — 

16 

UJ 

Print results 


[ Reinitialization 

i 


M 

r 


End 


Figure 2L - Flow chart for monopropellant system. 


Computer Program 

The formulation of the computer program follows the technique outlined in refer- 
ence 8. A flow chart of the program for linear compliance at the pump and injector is 
given in figure 21. Very minor changes are necessary if the flow chart is used for hyper- 
bolic accumulations. An explanation of the steps shown in the flow chart in figure 21 is as 
follows; each of the following steps is correspondingly numbered in the flow chart: 

(1) Read in the following initial conditions and physical constants as data: 

LI number of tank outlet component 

L2 number of pump component 

L3 number of injector component 

K1 number of time intervals for wave travel between friction orifices in tank 

K2 number of time intervals for wave travel between friction orifices in feed 

line 
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K3 number of time intervals for wave travel between friction orifices in 

discharge line 

K4 largest of K terms 

K5 number of time intervals for combustion dead time 

CT sonic velocity in tank, ft/sec 

CP sonic velocity in feed line, ft/sec 

CD sonic velocity in discharge line, ft/sec 

VAA amplitude of sinusoidal velocity perturbation of tank inlet, ft/sec 

VBB amplitude of sinusoidal velocity perturbation of pump, ft/sec 

R1 ratio of feed line to tank area 

R2 ratio of feed line to discharge line area 

2 

AF area feed line, ft 

A, B, C static pump curve coefficients 
BI coefficient of injector, ft^Vsec 

AT working time interval, sec 

M number of time intervals computations are to be carried out 

S frequency, cps 

BF1 friction orifice coefficient in tank, ft^^/sec 

BF2 friction orifice coefficient in feed line, ft 1 ^ 2 /sec 

BF3 friction orifice coefficient in discharge line, ft 2 /sec 

Q mean-line discharge, ft /sec 

PR1 ullage pressure, ft 

DH1 pressure change between tank friction orifices, ft 

DH2 pressure change between feed line friction orifices, ft 

DH3 pressure change between discharge line friction orifices, ft 

AP pump inlet compliance/feed line area 

APP injector dome compliance/discharge line area 

x reciprocal of gas residence time, 1/sec 

2 

y chamber constant, 1/ft 
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(2) Compute initial conditions and make assignments: 

Initial Conditions Assignments 

V1(L2,0) L=2, 1,L3 L = 1, 1, L3 

V0(L2, 0) HL(L, 0) J = -K4, 1, 0 

VI< L3 , 0) HR(L, 0) DHR(L, J) = 0 

V0(L3, 0) V(L,0) (L 4 L2,L3) DHL(L, J) = 0 

(3) Set up time loop and carry out computations for M cycles. 

(4) Assign a negative reflection for all waves reflected from liquid surface, 

(5) Set up loop for tank friction. 

(6) Set conditions and enter friction subroutine for tank. 

(7) Set conditions and enter subroutine tank. 

(8) Set up loop for feed line friction. 

(9) Set conditions and enter friction subroutine for feed line. 

( 10) Set conditions and enter subroutine pump. 

(11) Set up loop for discharge line friction. 

(12) Set conditions and enter friction subroutine for discharge line. 

(13) Set conditions and enter subroutine injector. 

(14) End time loop. 

(15) Print values of pressures and flows throughout the system. 

(16) Reinitialize as desired after a minimum of K4 intervals. 

The following numerical values for the computer program data were used in the calcu- 


lations: 


LI = 2 

CT = 1000 ft/sec 

L2 = 7 

CF = 2000 ft/sec 

L3 = 9 

CD = 3000 ft/sec 

K1 = 1 

R1 = 0.01 

K2 = 4 

R2 = 4 

K3 = 1 

AF = 0. 2 ft 2 

K4 = 4 

VAA = 0. 32 ft/sec 

K5 = 3 

VBB = 0.80 ft/sec 

A = 1208 

BI = 2. 8 

B = 668 

T = 0. 001 sec 
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C = -109 BFI 2 = negligible 

Q = 4 ft 3 /sec BF2 = 10 

PR1 = 172. 29 ft BF3 = 8 

DH1 = 0 

BH2 = 0 

DH3 = 0 

AP = 0, 0. 0001, 0. 00005, 0. 00001 
APP = 0,0.0001,0.00005,0.00001 
S= 5,7,9, 10,11, 12, 13, 14,16 

A listing of the FORTRAN IV program for linear compliance at the pump and injector is 
as follows: 


2 The tank friction was assumed to be negligible, therefore, no friction orifices were 
inserted in the tank. 
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COMMON/SOXl/AP.APP.AD.AF, t*x»y.at ,bo 

C0MM0N/80X2/C1 *C2 ,R 

COMMON /80X3 /A »B ,C *G»BF * ERROR 

G*32. 

DIMENSION PL (20*101) *PR ( 20*101 ) * V( 20.101 ) *VI ( 2*101 ) * 

1 VO (2*101) »DPR(20*133) »OPL(20.133) *VA( 101) *VT ( 101 > .00(20*133) » 

2YY2 ( 6000) 

C READ PROGRAM CONSTANTS 

1 READ (5*101) L1*L2*L3»K1*K2*K3*K4»K5,M»L5»L6»L7*L8 »KDDE2 
LL4*L2+1 
LL5=L3-1 
LLL=L1-1 
LL3=L2— 1 
KKK * 100 
K%2 s ICfCK + 1 

READ(S.IOO) T*S.BO, CT.CF.CD. VAA.VSB.Rl ,R2*AF,A,B,CC.8F1 ,BF 

12.8F3.DH1 .DH2.DH3. Q,AP ,APP , ERROR. X .DP 

100 FORMAT ( 8F10* 5 ) 

101 FORMAT (1515) 

600 FORMAT (5E16. 4) 

IF ( KQDF2 ) 2222.2221 *2222 
2222 READ( 5 . 100 )S 

IF(S.EQ..O) GO TO 1 
2221 LL2=L1+1 
MN * M 
KL2=-1 
AD*AF/R2 
AT=AF/Rl 

PR ( 1 . 1) *( 0/ AD/BO ) **2 - 225 *+9* *(Q/ (AD*BF3) )**2 

C INITIAL CONDITIONS IN TANK 2 TO Ll-1 

IFILLL-2) 90*91,91 
91 DO 60 L=2»LLL 
V(L.l ) =Q/AT 
PL(L*1)=PR(L-1 .1 )+DHl 

60 PR(L»1)*PL( L ,1)-(V(L,1 >/BFl)**2 

C INITIAL CONDITIONS AT TANK— LI 

90 PL(L1*1>=PR(LLL.1)+DH1 
VI ( 1 » 1 ) «Q/AT 
V0(1,1 J-Q/AF 

PRfLl.l )»PL(L1.1 )+(VK 1,1)*VI( 1,1)-V0( 1.1)*V0( 1.1))/(2**G) 

C INITIAL CONDITIONS IN FEEDLINE- Ll+1 TO L2-1 

DO 61 L*LL2»LL3 

V(L *1) =Q/AF 

PL ( L *1 ) “PR ( L-l » 1 1+0H2 

61 PR ( L » 1 ) *PL ( L ♦ 1 ) - ( V(L*1)/BF2) **2 

C INITIAL CONDITIONS AT PUMP- L2 

L*L2 ‘ 

PL ( L» l ) *PR (L— 1 *1 I+DH2 
VI ( 2* 1 ) =Q/AF 
VO ( 2*1 ) ®Q/AD 

PR( L» 1 ) S PL (L ,1 ) +A+B*0+CC*0*Q 
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C INITIAL CONDITIONS IN DISCHARGE LINE- L2+1 TO L3-1 

DO 62 L*LL4, LL5 
PL(L*1>»PR(L-1 »1 I+DH3 
V(L»1 )=Q/AD 

62 PR ( L * 1 ) =PL ( L *1 )-(V(L»l ) /BF3 )**2 

C INITIAL CONDITIONS AT INJECTOR- L3 

L*L3. 

PL ( L * 1 ) “PR (L-l *1 )+DH3 
DO 63 J*1 * K5 

63 QQILtJI'Q 
V(L3,1)*Q/AD 

PR I L » 1 1 = PL (L * 1)~(V(L»1)/B0)**2 
PEL=PR ( L » 1 ) 

Y=X#PEL/0 

c WRITE PROGRAM CONSTANTS AND INITIAL CONDITIONS 

WRITE! 6 *115) 

115 FORMAT I'1H1*6X*1HT*11X»1HS*1 IX »2H80, 10X »2HCT » 10X * 2HCP ,10X *2HCD. 10X 
WRITER *110) T »S*B0*CT »CF »CD» VA A »VB8 

WR I T£ ( £ * 116 > 

116 FORMAT ( 1HO, 6X, 2HAT, 10X.2HAF, 10X. 2HAD, 10X.1HA, 11X »1HB .1 IX »1HC,UX* 

1 3H8F1 «9Xi 3HBF2 , 9X *3H8F3 »9X»3HDH11 

13HBF1 ,9X»3HBF2.9X.3HBF3,9X.3HDH1 ) 

WRITE (6*117) 

WRITE (6 *110) 0H2 * DH3 *Q * AP ,APP * ERROR 

117 FORMAT {1HO»6X*3HDH2.9X»3HDH3«9X*1HO*1 IX *2HAP* 10X*3HAPP»9X*5HERROR) 
WRITE (6*110) AT »AF .AD.A,S*CC»8F1 ,BF2 »BF3 *DH1 

WRI TE<6*111 )L1*L2*L3*X1»K2*X3 *K4»K5»M,L5,L6 *L7.L8 
111 FORMAT ( 4H0L1® » I 5 * 5X * 3HL2= * I 5* 5X . 3HL3* * I 5 »5X, 3HK1 = * I 5 *5X *3HX2= * I 5 . 
15X*3HK3*. I5*5X,3HX4=»I5»5X*3HX5=» I5.5X,2HM=, I 5/ 4H0L5*, I5,5X»3HL6 

1**J5.5X*3HL7**I5.5X,3HL8=»I5) 

WRI TE( 6*502 ) 

WRITE (6 *501 KPL(L»1 ) .PR (L.l ) »V(L*1 ) *L *L = 1 »L3) 

502 FORMAT ( lHO »5X *2HPL »13X »2HPR , 1 3X . 1 HV * 1 OX . 1HL ) 

501 FORMAT (3F15. 5* !5) 

WRITE (6* 503) V0<1,1) ,V! (1 *1) *V0(2»1) ,VI(2*1) 

503 FORMAT ( AHOVO 3 *612,4,5X *3HVI * » G12 « 4 ) 

WRITE(6»120) 

120 FORMAT { 1H0»7X.5HPL-L1 . AX * 5HPR-L 1 * 4X*5HPL-L2» AX * 5HPR-L2 »4X *5HPL-L3 . 
1 4X*5HPR-L3 *4X* 5HVI-L2 ,4X *5HVI-L1 *4X *4HV-L3 • 5X * 5HPR-L5 »4X » 5HPR-L 
26, 4X»5HPR-L7*4X»2HB0,4X»3HVGP> 


FORTRAN ACTUAL MEANING 

PL ( L .1) * PL (L ,0) 

PR ( L ,1) = PR(L ,0) 

DPR(L»1 >*DPR(L.~K4) 

DPL (L»1)«DPL(L» ”K4 ) 
V(L»1)*V(L*0) 

VI ( 1 » 1 ) * VI (LI *0) 

VO (1*1) * VO ( LI » 0 ) 

VI(2»1) * VKL2.0) 

VO ( 2 » 1 ) * VO (L2 ,0) 

DO 2 L* 1 »L3 
DO 2 J=1»K4 

DPL ( L * J ) = .0 
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2 OPR (L*J) = .0 

CON = 2.*3.141593*S*T 

JK“K4— KI 

JL=K4~K2 

jM=K4~K3 

JN=K4-K5 

XX1*PR ( LX » 1 ) 

X2*PL ( L 1 * 1 ) 

X3 = 00 ( L3 * 1 ) / AD 

X4=PRtL5,ll 

XS=PL(L2*1) 

X6=PR( L2*l I 
XTsV! (1*1) 

X9=V ( L3 * 1 ) 

XlO=PL <L3»1 ) 

X11=PR!L2.1I 
Xl2=PR t L3 1 1 S 
X13=PR ( L5 * 1 ) 

X14=PR(L6»1) 

X 1 5 =PP f L7 * 1) 

X16=PR(L8. 1) 

PR (L3tl J=PR(L3 .1 )+DP 

C START TIME LOOP 

19 DO 3 J= 1 *MN 
J J= J+K4 
JJ J= J+K 5 
AJ = J+KL2+1 
JA= AJ 
J3 = JM-FJ 
J5= JN+ J 
J1=JK+J 
J2= JL+ J 

C REFLECTION FROM FREE LIQUID SURFACE 

DPR 1 1 ,JJ) =-DPL(2»Jl ) 

C FRICTION ORFICE ROUTINE FOR TANK RESISTANCE 

C=CT 

SF=BF1 

1FILL-2) 92,93*93 
93 DO 6 L=2,LLL 
112 FORMAT! 5F16. 4) 

6 CALL FOR(PLIL»J) *DPR<L-1* J1 I ,DPLIL+1, J1 ) »PR(L »J) ,VtL»J) ,PL IL*J+lI 
1 ,DPR (L*JJI*DPL ( L * J J I *PR (L*J+1I»V IL*J+1)I 
92 VT(J) a VAA*$IN(C0N*AJI 

C TANK ROUTINE 

L=L1 
Cl = CT 
C2 = CF 
R=R1 

CALL TANK (VI (1*J)*V0(1,J)*PL(L,J)»PR(L,JI *DPR (L-l , J1 ) ,DPL(L+1* J2 } 
l.VI (1»J+1 ) * VO (1*J+1) *PL(L,J+1I ,PR(L*J+1 ) ,DPL ( L , JJ > ,DPR(L,JJ) * VT ( J ) ) 
2) 

C FRICTION ORFICE ROUTINE FOR FEEDLIN6 RESISTANCE 
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C=CF 

BF=BF2 

DO 7L=LL2*LL3 

7 CALL F0R(PL(L*J) »DPR(L-1*J2I* DPL ( L+l ♦ J2 ) «PR (L»J)*V(L*J) *PL (L*J+1> 
1 »DPR (L»JJ),DPL ( L » JJ ) *PR (L.J+1>»V (L.J+H) 

C PUMP ROUTINE 

VA { J J = VB8*S1N(CDN*AJJ 

L S L2 

R=R2 

C1*CF 

C*CC 

C2-CD 

CALL PUMP (VA( J) »VI (2»J> *VOf 2* J> *PLIL»J) *PR(L»J> »DPR f L-l » J2 ) *DPL < L+ 
11 *J3» »Vl I2.J+1 ) * VO ( 2 * J+l > »PL { L * J+l ) *PR(L,J+1) ,DPLtL.JJ) » DPR<L*JJ> 
2) 

C FRICTION ORFICE ROUTINE FOR DISCHARGE LINE RESISTANCE 

C=CD 

BF=BF3 

DO B L = LL4 ,LL5 

8 CALL FOR I PL < L * J ) *DPR<L-1*J3> * DPL I L+l » J3 ) »PR( L*J) *V( L * J) *PL <L*J+1> 
l.DPR ( L * JJ } »DPL ( L * J J ) *PR (L»J+1)*V <L*J+1M 

C INJECTOR ROUTINE 

L=L3 

CALL I N JECT ( PL ( L » J ) »PR t L » J ) . V < L . J ) »QQ < L ♦ J ) ,DPR (L-l , J3 ) ,PL(L»J+1>* 
1 PR ( L » J+l ) *VIL *J+1 ) »QQ f L * JJJ ) ♦DPL(L»JJ) I 

C WRITE SOLUTION 

Y10=PL (L3.J+1 > -XI 
Y12=PR {L3.J+1 )-Xl2 
Y 1 3 = PR t L5 » J+l I-X13 
Y 1 4=PR { L6 »J+1 5-X1A 
Y1 S=PR { L7 * J+ 1 ) —X 1 5 
Y16=PR (LB »J+1 J-X16 
JA= AJ 

Y 1 = PR { L 1 * J + l ) - XXI 
Y2-PLIL1 » J + l ) - X? 

Y5=PL(L2»J+1)-X5 
Y6-PR(L2. J+l t-X6 
Y3=QQ<L3» J+l ) / AD-X3 
Y7=VI I 1 *J+1»-X7 
Y9=V(L3«J+lt-X9 
YY2 ( JA I =Y1 2 

WRITE <6* 108 ) JA.Y2 . Y 1 , Y5 . Y6 » Y1 Q , Y1 2 *Y3 »Y7 *Y9. Y13 .VIA, Y15 
108 FORMAT (15* 6F9. 3 * 3F9.4 * 5F8 . 3 I 

I F C ERROR. GE. 100. I GO TO 17 
IF (J-KKK) 3*10*1 

C END TIME LOOP 

3 CONTINUE 
110 FORMAT! 10G1 3 *4 ) 

10 KL2 =XL2+KKX 
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c 


RE-INIALIZATION 


DO 15 L«1 »L3 

V1L.11 ■ V1L.KX2) 

DO 65 J=1 .K5 
XXX1*XXX+J 

65 QQ( L3» J > S QQ{ L3 »XKK1 > 

DO 15 J=1 »X4 

XX = XXX+J 

DPL fl.J} = DPLIL.XXl 

15 DPR ( L » J 1 * DPR (L.XX } 

DO 16 L=1»L3 

PL < L » 1 1 * PL(L*XX2) 

16 PR <L.l) = PR < L »XX2 ) 
vm.ii * vki.xx2) 

V1I2.1I = VH2.XK2) 
VO(l.i) = VO ( 1 »XX2 ) 

VO ( 2 .1 ) * VO ( 2 »KK2 ) 

MN * MM-XXX 
IF(MN) 17.17.19 

17 IF< XODE2.NE.O) GO TO 2222 
60 TO 1 

END 
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